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^ ABSTRACT 

^0 We present the first numerical simulations of the solar interior to exhibit a 



tachocline consistent with the Gough & Mclntyre (1998) model. We find non- 
linear, axisymmetric, steady-state numerical solutions in which: (1) a large- 
scale primordial field is confined within the radiation zone by downwelling 
Cd meridional flows that are gyroscopically pumped in the convection zone (2) 



the radiation zone is in almost-uniform rotation, with a rotation rate con- 

^ sistent with observations (3) the bulk of the tachocline is magnetic free, in 

t^ 

^^^ thermal-wind balance and in thermal equilibrium and (4) the interaction be- 

^^ 

^f\ tween the field and the flows takes place within a very thin magnetic boundary 

■^ layer, the tachopause, located at the bottom of the tachocline. We show that 

O 

f^ the thickness of the tachocline scales with the amplitude of the meridional 

. . flows exactly as predicted by Gough & Mclntyre. We also determine the pa- 

> 

•^ rameter conditions under which such solutions can be obtained, and provide 

5^ a simple explanation for the failure of previous numerical attempts at repro- 



ducing the Gough & Mclntyre model. Finally, we discuss the implications of 
our findings for future numerical models of the solar interior, and for future 
observations of the Sun and other stars. 
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1 INTRODUCTION 

The rotation profile of tlie solar interior, which is nearly uniform within the radiation zone 



yet strongly differential within the convection zone ( Christensen-Dalsgaard & Schou 1988 

* E-mail: pgaraud@ucsc.edu (PG) 



2 L. A. Acevedo-Arreguin, P. Garaud & T. S. Wood 



Brown et al. 1989), presents a serious challenge to theoreticians (see reviews by Zahn 2007 



Garaud 2007). In particular, the thinness of the tachocline between these two zones implies 



that the transport of angular momentum below the radiative-convective interface must be 
predominantly latitudinal, yet the absence of rapid core rotation also requires significant 
radial transport throughout the radiation zone. At a cursory glance, this appears to be at 
odds with helioseismic sound-speed inversions, which suggest that the bulk of the radiation 
zone undergoes very little compositional mixing ( Elliott &: Gough]|1999 ). This apparent con- 
tradiction can be resolved if the radiation zone harbors a global-scale magnetic field, which 
transports angular momentum without engendering significant mixing. Such a field can also 



explain the uniform rotation of the radiation zone (Ferraro 1937 Mestel & Weiss 1987 



Riidiger & Kitchatinov 1997), as long as it remains confined strictly below the tachocline 



(Gough & Mclntyre 1998 MacGregor & Charbonneau 1999). The key problem to develop- 



ing a self-consistent model of the solar interior within this framework is to explain how the 
confinement of this field is maintained. 



The first model to address the "magnetic confinement problem" was proposed by Gough 



& Mclntyre (1998, GM98 hereafter). They argued that the field is confined by a large-scale 
meridional circulation, gyroscopically pumpecQby the convection zone, which converges at 
high latitudes and downwells into the radiation zone. In this model, the downwelling flows 
confine the magnetic field across a very thin layer located, by construction, at the bottom 
of the tachocline, which they called the "tachopause" . The thickness of the tachopause is 
determined by a balance between the downward burrowing of the flows and the upward 
diffusion of the field. At the same time, the remaining shear in the tachopause winds up the 
confined field lines, providing a Lorentz torque that pumps the meridional flow equatorward. 
The flows returns to the convection zone in a mid-latitude upwelling region, whose dynamics 
were not explicitly addressed in the GM98 model. 

In this picture, the bulk of the tachocline (excluding the tachopause) is "magnetic free" . 
Its dynamics are regulated by a balance of forces between thermal buoyancy, pressure, and 



the Coriolis force, leading to the so-called thermal- wind relation (Pedlosky 1979). The dif 



ferential rotation of the tachocline thus implies the presence of latitudinal temperature and 



entropy gradients, with a "hot spot" at the pole (GM98; Rempel 2005 Miesch et al. 2006) 



^ Gyroscopic pumping is a mechanism wliereby any azimuthal forcing, by conservation of angular momentum, also drives 
meridional flows. Gyroscopic pumping occurs in the solar interior by way of the turbulent stresses exerted by the convection 
zone onto the radiation zone. 
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The diffusion of the polar temperature anomaly is balanced by the advection of the back- 
ground stratification by the downwelling flow. In this way, the stratification and thickness 
of the tachocline constrain the mass flux allowed through it. Since the same mass flux must 
flow through the tachopause, and is regulated by the Lorentz torque therein (see above 
discussion), the thickness of the tachocline is tied to the strength of the internal magnetic 
field. 

Since 1998, a number of attempts have been made to obtain the GM98 tachocline in a 
self-consistent numerical model. Among them are the time- dependent, either axisymmetric 



or fully three-dimensional (3D) global numerical simulations of Brun & Zahn (2006), Rogers 



(2011) and Strugarek et al. (2011). These simulations are initialized with a global-scale 
dipolar magnetic field confined within a uniformly rotating radiation zone, and follow the 
resulting interaction of the field with differential rotation and meridional flows driven in the 



overlying layers. Brun & Zahn (2006) used a 3D model, but their computational domain only 
included the radiation zone, the top of which was modeled as an impenetrable boundary with 



an imposed latitudinal differential rotation. The more recent model of Rogers (2011 ) includes 



both the radiation and convection zones, but in a 2D meridional slice. Lastly, [Strugarek et al. 



(2011) have extended the 3D work of Brun & Zahn (2006) by including the convection zone 
in their computational domain. 

In all cases, the initially confined magnetic field gradually connects to the convection 
zone by magnetic diffusion. As a result, none of these studies recovers the picture of the 



solar interior envisaged by GM98. In the models of Brun & Zahn (2006) and Strugarek et al. 



(2011 ), the differential rotation propagates rapidly into the radiation zone once the field lines 



connect to the convection zone, as expected from Ferraro's isorotation law (Ferraro 1937). 



The resulting angular- velocity profile then bears little resemblance with observations. Rogers 



(2011 ) also finds that the field ultimately spreads throughout the whole domain, but that the 
radiative region remains mostly in solid-body rotation. The absence of Ferraro's isorotation 
probably results from the fact that the magnetic field is relatively weak in that model. Indeed, 
a relevant measure of the field strength is the Elsasser number, defined as A = B'^/AirprjQ 
(where B is the field amplitude, p is the density, rj is the magnetic diffusivity and Q is the 



rotation rate). In the simulations of Rogers (2011), A < 1 throughout the radiation zone, 
which may explain why the magnetic field appears to play no dynamical role. 

The failure of global numerical simulations to obtain solar-like dynamics have led some 
to conclude that a primordial magnetic field cannot explain the uniform rotation of the solar 



4 L. A. Acevedo-Arreguin, P. Garaud & T. S. Wood 



interior (Strugarek et al. 2011). However, we believe that this conclusion is premature. The 
numerical model parameters used in all simulations to date are necessarily very far from their 
corresponding solar values, because of computational limitations. In particular, the values 
used for the magnetic diffusivity and viscosity are typically many orders of magnitude larger 
than their microscopic solar counterparts. In each of the simulations reported above, the 
transport of magnetic field and angular momentum is dominated by these (artificially large) 
diffusivities. In the original model of GM98, on the other hand, viscosity is assumed to be 
entirely negligible, and magnetic diffusion is important only within the very thin tachopause. 
One should therefore consider very carefully the physical parameter conditions under which 
the GM98 model might be expected to apply. As we now show, this depends not only on 
the absolute magnitude of the diffusivities, but also on their ratios. 

As demonstrated by Wood & Brummell (2012[) (see also Garaud & Brummell 2008 



Garaud & Acevedo-Arreguin 2009 Garaud & Garaud 2008), the importance of viscosity in 



the radiation zone can be described in terms of the dimensionless parameter a, which is 
defined as 



o 



V N 



(1) 



where z/ is the viscosity, k, is the thermal diffusivity, N is the Brunt-Vaisala frequency, 
and ^0 is the Sun's mean rotation rate. This parameter can be interpreted as the ratio of 
the timescales for angular-momentum transport by Eddington-Sweet meridional flows and 
by viscous stresses respectively, across the same region. An analogous parameter appears 



in geophysical studies of stratified rotating flows (e.g. Lineykin 1955 Barcilon & Pedlosky 



1967). 



The condition for viscous effects to be negligible in a magnetic-free tachocline of thickness 
A is A ^ rcz/cT, where Tcz is the radius of the base of the convection zone. It is satisfied in the 



solar interior, where a varies from to about 0.5 within the tachocline (Garaud & Acevedo- 



Arreguin||2009), but is demonstrably not satisfied in any of the global numerical simulations 



described above. Indeed, the latter use realistic values for Qq and N, but have unrealistic 
diffusivities (and more importantly, an unrealistically large Prandtl number Pr = ly/n), 
leading to values of a that are significantly larger than one. We demonstrate in this paper 
that, in the regime a ^ 1, meridional flows downwelling from the convection zone are 
strongly suppressed and therefore unable to confine an interior magnetic field. 

There are essentially two different routes toward achieving the a < 1 regime in simula- 
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tions: by decreasing Pr while keeping N/Qq constant, or by decreasing N/^q while keeping 
Pr constant. The studies described above all took the first route, using solar profiles for 
iV/f^Q and making Pr as small as computationally possible. Unfortunately, to achieve a < 1 
in the tachocline would then require Pr < 10~^, which is far beyond the capabilities of any 
2.5D or 3D numerical code. Here, we favor the second route: we make Pr as small as possible 
in our numerical model and artificially decrease N/^q to achieve a < 1. Although in this 
approach neither N/^q nor Pr take their "true" solar values, we argue that this shortcoming 
is superseded by the need to achieve a non- viscous dynamical regime. 

A "proof-of-concept" of this approach was recently presented by IWood & Brummell 



(2012), using 3D numerical simulations in a local Cartesian domain straddling the radiative- 
convective interface. They show that an angular-velocity shear forced in the convection zone 
propagates into the radiation zone either by viscous diffusion or by advection by meridional 
flows, depending on whether a > 1 or a < 1. In cases with a > 1, they find that meridional 
flows decay exponentially with depth beneath the radiative-convective interface on a length- 



scale ~ Tcz/cT, as predicted by Garaud & Brummell (2008) and Garaud & Acevedo-Arreguin 



(2009), and that viscous stresses dominate angular-momentum transport. This explains why 
the global numerical simulations described above, which have a ^ 1 in the tachocline, are 
all dominated by viscous effects, and most likely also explains why the meridional flows 
are unable to confine the magnetic field. A similar conclusion had earlier been reached by 



Garaud & Garaud (2008 GG08 hereafter) albeit through more idealized simulations. By 



contrast, in cases with a < 1, Wood & Brummell (2012) find that viscous stresses are es- 
sentially negligible. Angular momentum is transported by large-scale meridional flows that 
burrow into the radiation zone on a local Eddington-Sweet timescale, as first discussed by 



Spiegel & Zahn (1992). These flows retain a significant amplitude in the radiation zone, as 



expected from the work of Garaud & Brummell (2008) and Garaud & Acevedo-Arreguin 



(2009). 



In short. Wood & Brummell (2012) showed that it is possible to obtain solar-like dynam- 
ics without using true solar parameters, by identifying the appropriate parameter regime - 
in the case of the tachocline this requires having a < 1. However, the numerical challenge 
involved in modeling simultaneously the convection zone and the radiation zone, in a full 
3D MHD simulation and in a parameter regime where a < 1, remains quite formidable. 
To prepare for a time in the not-so-distant future where such simulations are possible, we 
have already begun to study the problem using reduced models, which enable us to refine 
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our understanding of the dynamics of the system and provide insight into the appropriate 



selection of other parameters. In Wood et ah (2011), we studied a Cartesian and laminar 



toy model of the solar interior. Its analytical simplicity enabled us to understand in greater 
detail the structure and global force balance of the tachocline and the tachopause, leading 
us to propose a set of diffusion coefficients (viscosity, magnetic, and thermal diffusivity) 
that can realistically be used in simulations to yield a GM98-like tachocline. The next step 
is to use these proposed parameters in two concurrent and complementary models: in an 
axisymmetric, steady-state, nonlinear model of the full solar interior (this paper), and in a 
3D, fully nonlinear local Cartesian model of the tachocline (Wood & Brummell, in prep). 

In what follows, we consider a global quasi-steady axisymmetric numerical model of 
the solar interior, presented in detail in Section [2j Section [3] discusses an important yet 
subtle issue discovered in the process of searching for tachocline-like solutions, namely the 
difference between magnetic confinement of the type considered by GM98, which occurs at 
the bottom of the tachocline, and confinement by fiows in the convection zone, which occurs 
at the radiative-convective interface. This second type of magnetic confinement is similar 



to that proposed by Kitchatinov & Riidiger (2006), and depends sensitively on the details 



of the model at the bottom of the convection zone. In Section |4| we present and discuss the 
first numerical simulation of the solar interior to exhibit the layered tachocline/tachopause 
structure anticipated by GM98. Section |5] then examines how various properties of the 
tachocline and tachopause vary with magnetic and thermal diffusivities, and ultimately 
validates the predictions of the GM98 model. We also run a set of simulations in a parameter 



regime similar to that used by Strugarek et al. (2011), and recover their results - namely, 



that the magnetic field is not confined by the tachocline fiows - hence showing that using 
the correct value of a is indeed necessary if one wishes to obtain a GM98-like solution. 
Finally, we discuss the implications of our results in view of future 3D simulations of the 
solar interior, and within the greater context of stellar astrophysics in general, in Section |6| 



2 THE MODEL 

Our goal is to study how the upper layers of the radiation zone respond to forcing by 
the differential rotation in the convection zone above and to the presence of a large-scale 
primordial field below. To do so, we use a model and numerical algorithm that is an extension 
of the one developed by GG08. 
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We focus primarily on getting the dynamics of the radiative interior (radiation zone and 
tachochne) right, at the cost of having to simphfy the dynamics of the convection zone quite 
dramatically. This sacrifice is necessary, as no numerical algorithm today can model at the 
same time the very fast timescales and short lengthscales associated with convection, and 
the very long timescales and global lengthscales appropriate within the deep interior. Our 
convection zone model is described in detail in Section 12.5.11 

The timescale of principal interest in the magnetic confinement problem is the timescale 
for magnetic diffusion across the tachocline, which in the Sun is approximately 10 Myr. This 
is much shorter than the global-scale Eddington-Sweet and magnetic diffusion times, but 
much longer than the oscillation periods of inertia, gravity, and Alfven waves. On this 10 
Myr timescale, we may assume the tachocline to be in a quasi-steady force and thermal 
balance, as in GM98. The effect of the fast dynamics (including the mean effects of waves 
and turbulence) can be parameterized, if desired, through the inclusion of additional terms 
(e.g. Reynolds stresses, turbulent thermal diffusivity, turbulent magnetic diffusivity, etc.). 
In this paper we use such parameterizations only within the convection zone (see below). 
This is for simplicity, and because the effect of fast dynamics on the tachocline is neither 
particularly well understood, not well constrained. We discuss their potential impact on our 
results in Section |6l 

The thermodynamic structure of the solar interior is accurately known from helioseis- 
mology, and is not greatly affected by the flows in the radiation zone. We therefore linearize 
all the governing equations around a hydrostatic background state derived from the 1-D 



solar model of Christensen-Dalsgaard (1996). We then solve numerically the resulting set of 
MHD equations, assuming axisymmetric dynamics. 

GG08 designed a numerical algorithm to search for steady-state solutions of this system 
for various input parameters, such as the assumed viscosity, thermal diffusivity, magnetic 
diffusivity, etc. In this work, we use the same algorithm, but with a few notable modifications. 
We now describe the original method and the modifications in detail for completeness. 



2.1 The background state 

As in GG08, our background state assumes a spherically symmetric Sun in hydrostatic equi- 
librium. All quantities are expressed in spherical coordinates {r,6,(j)), where 6 is colatitude. 
Whereas GG08 modeled the radiation zone only, with the radial coordinate r ranging from 
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Tin = O.35-R0 to Tout = O.7I-R0, we now include a significant portion of the convection zone as 
well in our numerical domain, and probe deeper towards the center. In what follows, r spans 
the interval [rin = O.O5-R0,rout = O.Oi?©]. We interpolate Model S (Christensen-Dalsgaard 



1996) in that range to obtain the density p, temperature T, pressure p, gravity g, and heat 



capacity at constant pressure c„ of the background fluid as a function of r. 



2.2 The model equations 

We work in a frame that rotates with angular velocity 
a2 3a4 



1^0 = u 



eq 



35 



(2) 



about the vertical axis, where f^cq = 27r x 463nHz is the equatorial velocity of the Sun near 



the base of the convection zone, 02 = 0.17, and 04 = 0.08 (Schou et al. 1998 Gough 2007) 



In this frame, the total specific angular momentum of the convection zone is zero ( Gilman 



et al. 1989). 



We first expand each of the thermodynamical variables q as q{r,6) = q{r) + q{r,6), 
where the bar indicates the spherically symmetric background state and the tilde refers to 
the axisymmetric perturbation. Under the steady-state assumption and once linearized in 
the thermodynamical fields around the background state, the system of governing equations 
can be written as: 

p(r) 



2p{r)Q,Q X u + p(r)u ■ Vu 
p(r)T(r)u ■ Vs(r) 

V X (u X B) 
P 



-Vp - pg{r)er + j X B + V ■ n 

V ■ {k{r)Vf) 

Vx [r/(r)(V xB-47rjo)] 



(r) 



U — Ur^ 



p{r) 
V ■ (p(r)u) 

VB 



+ 



T 



-p{r) T{r) 








(3) 
(4) 
(5) 
(6) 
(7) 
(8) 



where u is the velocity of the fluid, B the magnetic field, j the electric current density, 
7] the magnetic diffusivity, and k the thermal conductivity. The viscous stress tensor 11 
incorporates the contribution of the viscosity v. 



n = p{r)v{r) 



Vu+(Vu)^--(V-u)I 



(9) 



where I is the identity matrix. Note that the advection of the entropy perturbations has 
been neglected. Finally, and as in GG08, we retain only those terms in u ■ Vu that include 
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the azimuthal velocity u^, since the other (meridional) velocity components are generally 
much weaker. We do this to accelerate numerical convergence, but have checked that this 
has very little effect on the final solution in the parameter regime of interest. 

This system of equations is very similar to the one used by GG08, with two notable 
differences. First, we parametrize the driving of the differential rotation in the convection 
zone by the body-force term — -(u — Ucz) where 1/T{r) vanishes in the radiation zone (see 



Section 2.5.1 for more detail). Second, we have introduced a source term 4'7rjo = V x Bq 



in the magnetic induction equation to maintain an assumed primordial magnetic field Bq 



against diffusion (see Section 2.5.2 for detail). 



2.3 The diffusivity profiles 

Although our steady-state model is able to achieve much lower values for the various diffu- 
sivities (viscosity z/, magnetic diffusivity rj and thermal diffusivity k, = k/pcp) than direct 
numerical simulations, they must still be artificially increased by several orders of magnitude 
to yield a numerically tractable system. As a result, there is little benefit in using realistic 
profiles for these quantities, so we choose them instead to be as simple as possible. 

In what follows, we use Tcz = O.7127-R0 to denote the radius of the radiative-convective 
interface. For simplicity, we take the diffusivities to be constant within the radiation zone 
(r < Tcz), with values Uj.^, rjj-z, and /tj-z respectively. In the convection zone, we model the 
effects of the turbulence on the magnetic field and heat transport as an increase in the 
effective diffusivitiesjj For reasons that are explained in Section |3| we construct the profiles 
to ensure that this increase occurs slightly above the radiative-convective interface, at a 
radius ri > rcz. The choice of ri is discussed in Section |3]as well. Thus, 

u{r) = V,, (10) 

riir) = r^rz + H{r - ri)r]t{r) , (11) 

K(r) = Krz + H{r - ri)/tt(r) , (12) 



where H is a Heaviside function, and 

' r — r2 



Vt{r) = ^[Vcz-rir 



rout - n 



1 + tanh 



Ao 



(13) 



^ We do not increase the viscosity in the convection zone because we have already introduced a body-force to parameterize 
the turbulent transport of momentum. 
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Figure 1. Non-dimensional diffusivity profiles v/B3^Q,q (dashed line), k/Rq^q (dotted line) and tj/R^Qq (solid line) as a 
function of r for our reference model (see Section HI. The constant non-dimensional values in the radiation zone are Ei, = 
5.0 X 10-1", Er, = 1.5 X 10-9, and E^ = 5.0 x 10-^ 



Utir) 



I rLp7 Ihr 



r — n 

'^out - n 



1 + tanh 



r — r2 



(14) 



The tanh terms in (13) and (14) smooth the transition between the laminar and turbulent 



regions over a layer of thickness A2 centered at radius r2. In all that follows, we take r2 = 
Tcz -|- O.OSRq and A2 = O.O2_R0. The resulting profiles, as used in our reference model, are 
shown in Figure [TJ The selected values of ri, k^z and rjcz for this model are reported in Table 
1. While ri has significant influence on the solution (see Section |3] for detail), Kcz and r/cz do 
not as long as they are large enough. 



2.4 The parameter a and the selection of A^^ 

In Section hi we introduced and discussed the importance of the parameter a = \^FiN/Qq 
in setting the vertical scale of penetration of the meridional flows into the radiative interior. 
We therefore have to ensure that a is of the same order as in the Sun. Since our Prandtl 
number is not solar, we choose to modify the background stratification profile N in the 
simulation in such a way as to have a{r) be the actual solar profile (TqIt), as computed from 



Model S ( Christensen-Dalsgaard 1996). To do so, we proceed as follows. 



We first compute the "true" background diffusivities z/0(r) and K,Q{r) in the radiation 



zone from the formulae given by Gough (2007). We then compute the solar stratification 
parameter 



Dynamics of the solar tachodine 11 



where NqIt) is interpolated from Model S. Once the numerical diffusivity profiles z/(r) and 



K(r) have been selected (see Section 2.3), we construct an artificial background buoyancy 
frequency profile N satisfying: 

{ap, (r)ilp. . / — ^ for r < r^z 
°^ ^ °V^rz (16) 

otherwise. 

Figure ^ shows N{r) as used in our reference model, and compares it with the solar profile 
NqIt). In order to have a solar a, we have to take N{r) to be a factor of about 20 smaller 
than in the solar tachocline. 

As we demonstrate in Sections |4] and [5} using a realistic value of a is crucial to obtain 
confined-field solutions (see also Wood & Brummell 2013, in prep). But what are the conse- 
quences of choosing a lower-than-solar value of N on other aspects of the problem? In our 
steady-state, laminar model of the solar interior, choosing a lower value of N can only influ- 
ence the large-scale meridional flows, and does so in such a way as to improve the dynamical 
realism of the results. In 3D time-dependent simulations, however, a lower value of N would 
also increase the depth of penetration of overshooting convective plumes, and change the 
frequency of oscillation of gravity waves, which may have detrimental effects on the solu- 
tions. However, it is clear that no simulations can ever yield an exactly solar tachocline until 
such a time when they can be run at exactly solar parameters. In the meantime, one can 
still meaningfully study the dynamics of the tachocline even with a thicker overshoot layer, 
by selecting other input parameters to ensure that the tachocline is correspondingly thicker 
as well. 

2.5 The forcing terms 

The main differences between the model used here and that of GG08 lies in the forcing terms. 
In GG08, forcing was applied through boundary conditions. At the outer boundary of their 
domain, differential rotation and radial inflows/outflows were imposed to model the driving 
of meridional flows by differential rotation in the overlying convection zone. A primordial 
magnetic fleld was maintained by a point-dipole at r = 0. This method, however, has two 
disadvantages. First, it cannot take into account the back-reaction of the radiation zone on 
the convection zone dynamics, and introduces artiflcial Ekman layers at the outer boundary. 
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Figure 2. Background buoyancy frequency profile N{r) in our reference model (see Sectionl4jl, compared with the solar profile. 
A lower value of TV is needed to ensure that <j{r) is exactly solar in our numerical model. 

Secondly, under the point-dipole assumption, the primordial field Bq reaches unrealistically 
large amplitudes as r — )■ 0, leading to numerical convergence difficulties and driving artificial 
MHD instabilities near the lower boundary. Both problems ultimately prevented GG08 from 
obtaining solutions at low-enough diffusivities to exhibit a satisfactory GM98-like tachocline. 
In order to model the dynamics of the complete system, including both the convec- 
tion zone and the radiation zone, we have extended the computational domain to r G 
[0.05, O.9]i?0, as discussed earlier. In order to improve the stability of the numerical scheme, 
we have also modified GGOS's algorithm to drive the system through body-forcing terms 
rather than boundary conditions. These terms are now described in more detail. 



2. 5. 1 Convection zone forcing 

In the solar convection zone, gyroscopic pumping is caused by the Reynolds stresses associ- 
ated with strongly anisotropic eddies, which drive the system away from uniform rotation 



(Riidiger 1989; Kichatinov & Riidiger 1993 Rempel 2005). To model this in detail either 



requires the use of a 3D time-dependent code as in Strugarek et al. (2011 ) or Rogers (2011 ), 



where angular-momentum transport arises naturally from the convective dynamics, or the 
use of a parametric model for the convective Reynolds stresses (Kichatinov fc Riidiger||1993 



Rempel 2005; Garaud et al. 2010). Both approaches have well-known pros and cons. The 



former is computationally expensive, and of course cannot be implemented within the steady- 
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state axisynimetric approach we have elected to follow here. The latter, by contrast, can. 



and has been used with some degree of success by Rempel (2005) to model the differen- 
tial rotation profile and meridional flows within the convection zone, but its reliability is 
inherently tied to the reliability of the turbulence model used, a factor which is difficult to 
evaluate objectively. 

Here, we use a much simpler model for the Reynolds stresses, which takes the form of 



a Darcy friction law (Bretherton & Spiegel 1968 Garaud & Acevedo-Arreguin 2009 Wood 
er"aL][20lI| ): 



PJr) 
r{r) 



(u-Uc 



(17) 



where u^z is the observed azimuthal velocity of the solar convection zone (expressed in the 
rotating frame). 



u^ 



■ sin 9 [fieq(l — 02 cos^ 9 — a^ cos"^ 9) — QqJ e^ 



and where the relaxation timescale T(r) is defined such that 



1 



T{r) 



Tr. 



r — Tr 



'out ' c 



H(r — Vr 



(18) 



(19) 



The quantity Tcz, which we take to be O.lfig^, can be interpreted as a turnover time for the 
convective fiows. We have chosen the r profile so that 1/T{r) starts increasing linearly exactly 
from the base of the convection zone upward. By increasing the forcing continuously, we avoid 



the formation of artificial boundary layers at the radiative-convective interface ( Garaud & 



Acevedo-Arreguin 2009). In other respects, the particular choice of the function T{r) has 



relatively little impact on the resulting long-term dynamics in the radiative interior (Wood 



et al.|20lT), see below for detail. Finally, note that the forcing that maintains the differential 



rotation also produces a persistent gyroscopic pumping of meridional fiows in the convection 



zone ( [Garaud fc Acevedo-Arreguin||2009[ [Wood et al.||2011t [Wood fc Brummell||2012[ ), with 
downwelling in polar regions and upwelling near the equator. 



The reason for our choice of forcing is mere simplicity. However, as shown by Wood 



et al. (2011), the dynamics of the tachocline itself are fairly independent from those of 



the convective zone, as long as, to leading order, (1) the exchange of angular momentum 
between the convection and radiation zones is less efficient than the redistribution of angular 
momentum within the convection zone (2) the convective heat transport is efficient enough 
that the convection zone is very nearly adiabatic and (3) the meridional fiows within the 
tachocline remain significantly weaker than the meridional fiows in the convection zone. 
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These three conditions are naturally satisfied in the real convection zone, and were implicitly 



assumed in the models of Spiegel & Zahn (1992) and GM98. 



The first two conditions can be satisfied in our Darcy friction model provided Tcz^q < 1 
and Kcz > {Rq — fcz)^!^- This guides our choice of r, and Kcz, although their exact values 
have little-to-no infiuence on the resulting system dynamics as long as the inequalities above 
are satisfied. The third condition is automatically satisfied in the Sun, which has a relatively 
massive convection zone and strong differential rotation (both aspects are represented fairly 
accurately in our model), but may not be in higher-mass stars. 

2.5.2 Forcing of the primordial magnetic field 

The primordial field Bo is maintained against diffusion by a source term in the induction 
Equation (|5| that imposes a permanent azimuthal electric current jo = V x Bo/47r. To avoid 
artificially altering the tachocline dynamics, the imposed current is localized within a region 
'"a ^ ''" ^ ^b deep in the radiation zone: 

hir,e) = \ ^e '^ (20) 

I otherwise 

where r^ = O.IRq and rj, = O.SRq. In the absence of any induction by fiuid motions, this 

current would generate an open dipolar magnetic field Bq, whose exact expression is derived 

in Appendix |Aj We do not impose any current close to the core, in order to avoid the 

numerical difficulties described by GG08, and so Bq is uniform within the sphere r < r^. 

The amplitude of the current density, Jq, is chosen so that the magnitude of this uniform field 

is equal to Bq. Figure [3] shows selected magnetic field lines of Bo to illustrate its structure. 

Note that the actual structure of the poloidal field calculated from the simulation is only 

similar to Bo when the magnetic Reynolds number is much smaller than one. 

2.6 Boundary conditions 

The inner core (r < Tin, not modeled explicitly) is assumed to be a thermally and electrically 
conducting solid with the same thermal and magnetic diffusivity as the fluid just above 
r = Tin- As such, the bottom boundary at nn is impermeable. The solid core has a uniform 
angular velocity, i^jn, chosen in such a way as to guarantee that it exerts a zero total torque 
on the fluid above. Altogether, these conditions imply: 

Ur{rin,0) = (21) 
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0.0 



Figure 3. Geometry of the computational domain and of our assumed primordial magnetic field Bq. The latter is maintained 
by the electric current density jo (see Equation | |20[ l and Section |2.5.2| for detail). The solid curves represent a few selected field 
lines. The outermost circle (solid line) marks the solar surface. The dashed circles, from the outermost to the innermost one, 
mark rout, fez and ri^. The shaded area is the interval between ra and ri,, and marks the region where the azimuthal current 
jo is imposed. 



due{rin,0) 
dr 



tt/2 



-7r/2 



on 



BrB, 



r-D(j> 



pi'r^„ -TT- sin 9 H ri^ sin ^ | sin ^ d^ = . 

or 4:11 



(22) 
(23) 



The temperature perturbations in the core satisfy Laplace's equation: 

V^T = . (24) 



The solutions of Equation (24) are matched onto the solutions obtained in the computational 



domain at the boundary r = rin, requiring continuity of T and dT /dr. Similarly, the magnetic 
field in the core satisfies: 



V^B = 



(25) 



and the solutions in the computational domain are matched onto the core solution by re- 
quiring continuity of B^. and Bq at r = r^^. 

At the outer boundary, r = rout, we impose the following boundary conditions: 



Ur{r out, 0) = 

U<j,{r out, 6) = Mcz(^out,l 



(26) 
(27) 
(28) 
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which represent an impermeable, no-shp outer boundary, where Mcz is given in Equation ( 18 ). 
As with the inner boundary conditions, we assume that the temperature perturbations and 
the magnetic field satisfy Laplace's equation outside the domain, and require the continuity 
of T, dT jdr^ Br and Bg across the boundary. 



2.7 Non-dimensional parameters 

For comparison with other works, we now present a non-dimensional version of the governing 
parameters. By normalizing distances with respect to the solar radius Rq, time to the inverse 
of the mean solar rotation rate Qq, velocities to RqQq, density to po = Ig/cm^, and the 
amplitude of the magnetic field B to Bq (the constant value of the imposed primordial field 



in the core, see Section 2.5.2), the following non-dimensional parameters naturally arise: 



-C/K - 


RqHq 


Eri = 


Rq^Iq 


P — 


l^vz 


-C/l/ 


Rl^Q 


A = 


Bl 



(29) 

47rpo^rz^0 

where we take po = Ig/cm'^. The first three numbers are the "Ekman" numbers for the 
radiation zone, which are dimensionless measures of the three diffusivities. The last is an 
Elsasser number, which measures the relative strength of the Lorentz force and the Coriolis 
force associated with azimuthal fiows, in the core. The effective Elsasser number in the 
tachocline is typically significantly smaller than A, since the magnitude of Bq drops rapidly 
with r (see Appendix A for detail). 



2.8 The numerical method 

The system of partial differential equations and boundary conditions discussed in this Sec- 
tion is solved numerically using the nonlinear relaxation Newton-Raphson-Kantorovich al- 
gorithm developed by GG08. The details on how these equations are solved and implemented 
in our model are described in GG08 (see their Appendices). 
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3 MODEL SENSITIVITY TO TURBULENT MAGNETIC DIFFUSIVITY 
PROFILE IN THE CONVECTION ZONE 

In Sections |4] and [5} we study the model results in detail and discuss their implications for 
the dynamics of the solar tachocline and radiation zone. First, however, we must discuss an 
important point regarding the model sensitivity to the magnetic diffusivity profile near the 
radiative-convective interface. Figure |4] illustrates this issue by showing the results of two 
simulations that differ only in the location ri at which the magnetic diffusivitjjj increases 
from its laminar radiation zone value to its turbulent convection zone value (see Equation 



(13) in Section 2.3). On the left, ri = Tcz, while on the right, ri = r^z + O.OO2i?0. 

When ri = Tcz, we see that the field is unconfined at high latitudes, but when ri is 
moved up even a tiny distance into the convection zone then the magnetic field is confined 
by large-scale meridional flows at the radiative-convective interface. For the sake of brevity, 
in all that follows we call this effect "pre-confinement" , to contrast it with confinement by 
meridional flows deeper in the radiation zone. Pre-confinement is quite different from the 
dynamics described by GM98 for the reasons pointed out below. 

But first, in order to understand why such a small change in ri results in such a dramatic 
effect on the magnetic field configuration, we compute the magnetic Reynolds number Rm = 
\ur\L/ri, where L is (somewhat arbitrarily) selected to be L = O.IRq. Since the value of Ur 
in the lower convection zone is essentially identical in both cases (see Figure [sk), Rm is 



mostly controlled by the magnetic diffusivity profile. When ri = r^z, the flow velocity drops 



more-or-less at the same rate as rjt (see Equation (13)) while approaching the base of the 
convection zone, leading to Rm < 1 for all r > rcz, as shown in Figure |5)d. On the other 
hand, when ri —rcz = O.OO2i?0, a shallow region r^z <t < ri appears where the same flows, 
this time combined with a sufficiently low magnetic diffusivity, yield Rm ^ 1. Such flows 
can then substantially affect the magnetic field, advecting it horizontally just above the base 
of the convection zone, which results in the aforementioned magnetic "pre-confinement" . 
The possibility of magnetic confinement by latitudinal flows at the bottom of the con- 



vection zone has previously been discussed by Kitchatinov &; Riidiger (2006). They argued 



that viscosity (either microscopic or turbulent) prevents meridional flows from penetrating 



^ Changing ri also changes the position where the thermal diffusivity increases, see Equation jl4[ l; however, this has a negligible 
effect on this discussion, since the thickness of the thermal diffusion layer is much larger than the thickness of the magnetic 
diffusion layer. 
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Figure 4. Comparison of two simulations with and without prc-confinement (see Section Is] for detail), produced using Ei, = 
5 X 10~^, Erj = 1.5 X 10~*, and E^ = 5 x 10^'' (i.e. diffusivities that are all 10 times higher than the reference model studied 
in SectionWl, a solar a profile, and A = 0.66 X 10^. On the left ri = rcz, and on the right ri — rcz = 0.002ijQ. In both figures, 
starting from the top left panel and going clockwise, we show selected flow streamlines (with solid lines for clockwise flow, 
and dotted lines for anticlockwise flow), temperature perturbations (in units of Kelvin), angular velocity (in units of fieq), and 
magnetic field lines. The strip below each panel zooms into the region between r = 0.5Rq and rcz- 
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Figure 5. Left: Radial flow velocity \ur\ at 80° latitude, for the two simulations shown in FigureW] The solid line shows upwelling 
flows, and the patterned line shows downwelling flows (dotted for the ri = rcz case, and dashed for the ri — rcz = 0.002i?Q 
case). The vertical velocities in the convection zone are indistinguishable. Right: Magnetic Reynolds number, Rm = \ur\L/ri, 
computed at 80° latitude, where L = O.li?0. The same linestyle as the left-side flgure is used. The two kinks in the curves 
correspond to rcz = O.7127il0 and ri = O.7147R0. Just above the base of the convection zone, Rm 2> 1 for the case where 
''I > '"cz while Rm <SC 1 for case where ri = rcz- 
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more than a short distance into the radiation zone, leading to a strong radial gradient in 
the latitudinal flow velocity at the radiative-convective interface. This shear, combined with 
an assumed weak magnetic diffusivity, leads to field confinement much as in Figure lib. We 



note, however, that the model of Kitchatinov & Riidiger (2006) does not explicitly consider 



angular-momentum balance, instead treating magnetic confinement as a purely kinematic 
problem. It is therefore not a fully self-consistent model for the tachocline. Moreover, as 
discussed in Section [1} and illustrated in Figure |4| viscosity does not prevent meridional 
flows from penetrating into the solar radiation zone. 

What is interesting, however, is that pre-confinement clearly facilitates confinement 
deeper in the radiation zone, at least within the scope of our numerical model. With pre- 
confinement, we are able to find deeply confined solutions (see Figure lljo and Section 111), 
while similar solutions are much more elusive in models which are not pre-confinecQ 

Is magnetic field pre-confinement by the convection zone flows a possible scenario in the 
Sun? We believe it is, although not in the manner described above: in the solar convection 
zone, transport of magnetic field is probably dominated by the action of small-scale turbulent 
convective plumes, rather than advection by large-scale mean meridional flows. However, it is 



often argued that such turbulent transport also promotes magnetic confinement ( Zel'dovich 



1957 Radler 1968 Weiss 1966; Drobyshevski & Yuferev 1974). In the case of the solar 



tachocline, this has been demonstrated in a series of numerical simulations by Tobias et al. 



(2001) (see also Garaud & Rogers 2007) who showed that overshooting convective plumes 
produce a net pumping of magnetic field out of the convection zone. In other words, whether 
by large-scale laminar flows, or by small-scale turbulent flows, pre-confinement is likely to 
take place. 

Pre-confinement alone, however, cannot explain all tachocline-related observations. In- 
deed, if magnetic pumping were the only confinement mechanism in operation, then the 
tachocline would be only as deep as the convective overshoot layer, which is thought to 



be significantly smaller than the observed tachocline thickness ( Christensen-Dalsgaard et al 



1995 Brummell et al.|2002 Rogers et al.|2006 ); such a model is not consistent with helioseis- 
mology. Moreover, turbulent magnetic pumping is only known to be effective at confining 
a horizontal magnetic field. Recent simulations by Wood and Brummell (2013, in prep.) 



* This does not mean they don't exist. There are strong indications that they do at lower diffusivities than the ones for which 
we are able to obtain fully resolved, well-converged solutions. 
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Parameter Radiation zone (r < Tcz) Convection Zone (r > rcz) 



u/Rl^Q 


5 X 10-1° 


v/Rl^Q 


1.5 X 10-9 


k/RIUq 


5 X 10-'^ 


a 


o"©('') 


T 





A 


0.66 X lO** 



5 X 10-1° 

1.5 X 10-9 + 5.0^:^;^ [l + tanh (^)] H{r - n) 



5 X 10-^ + 0.5^:^;^ [l + tanh (^)l H{r - n) 


1 [" y — rcz 
Tea [rout-l-c! 

0.66 X 10« 



Table 1. Non-dimensional parameters for the reference model. In addition, ri — r^z = O.OOOS-Rq, and r2 and A2 were defined 
in Section |2.3| Finally, rczf2o =0.1. 

strongly suggest that large-scale flows similar to those proposed by GM98 are needed to 



confine the field near the poles, as previously proposed by Wood & Mclntyre (2011). 

In summary, since turbulent pre-confinement via magnetic pumping is likely to be present 
in the Sun, and since pre-confinement clearly facilitates our task of finding deeply confined 
solutions reminiscent of the GM98 model, in all that follows we let it take place (in the 
laminar sense) by choosing ri — r^z = O.OOOSi?© (an even smaller separation than that 
illustrated in Figures H^ and pi) unless otherwise indicated. 

4 REFERENCE MODEL RESULTS 

We now present and analyse in detail the results of our "reference" model. The governing 
parameters for this simulation were selected after a series of experiments starting from a 



parameter regime similar to the one proposed by Wood et al. (2011) and gradually reduc 



ing all diffusivities simultaneously. We found it necessary to reduce all diffusivities by four 



orders of magnitude, relative to those suggested by Wood et al. (2011), in order to obtain 
a solution with a well-formed tachocline and tachopause . The reference model parameters 
are summarized in Table [ij We first describe its qualitative properties and then analyze it 
more quantitatively by studying the balance of forces in the tachocline. 

4.1 A first glimpse of the solar tachocline 

Figure [6] shows the fiow field, temperature perturbation, magnetic field, and angular ve- 
locity in the reference model, demonstrating that magnetic confinement strictly below the 
radiative-convective interface, i.e. distinct from "pre-confinement" , is possible. 

The top left figure shows the steady-state meridional fiows in the reference model. From 
the surface to the center, we identify three main regions. In the convection zone, the flows 
are driven by the gyroscopic pumping effect of the imposed differential rotation. Fluid is 
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Figure 6. Steady-state solution for our reference model. The four panels show the same variables as in Figure 111 Large-scale 
meridional flows generated in the convection zone penetrate into the radiation zone at high latitudes. These downwelling 
flows are deflected by the internal magnetic field in the tachopause and, in turn, confine the field. The confined field enforces 
solid-body rotation of the radiative interior below the tachopause. 



pumped toward the poles near the surface, and returns equatorward near the base of the 
convection zone. The tachochne is located below the radiative-convective interface. Since 
the radiative region in our reference model is weakly stratified in terms of the parameter a, 
a fraction of the meridional mass flux driven in the convection zone enters the tachocline. 
There, we observe two cells. The main tachocline cell downwells at the poles, is deflected 
equatorward at about 0.5Rq, then returns to the convection zone in a thin upwelling region 
around 30° in latitude. Close to the equator, a small cell with meridional flows circulating 
in the opposite direction is also visible. The downwelling flows in the large cell, although 
weak in amplitude, are sufficiently strong to distort the magnetic field (see below). Below 
the tachocline, we observe a thin meridional counter-cell, which is part of the tachopause. 
Deeper within the radiation zone, meridional fiows become negligibly weak. 

Note that the tachocline in this simulation is much thicker than in the Sun. This is 
not surprising, since the model thermal and magnetic diffusivities are much larger than 
in the Sun, and we expect the thickness of the tachocline and the tachopause to depend 
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sensitively on these parameters. In Section [5| we test the dependence of this thickness on 
the diffusivities, and demonstrate that our results become consistent with observations once 
extrapolated to solar parameter values. 

The bottom left figure shows the effect of these flows on the poloidal magnetic field. In 
the upper part of the convection zone, the magnetic diffusivity is high and the flows do not 
affect the field. By contrast, we see that the field lines are strongly distorted by the flows in 
the lower part of the convection zone, where the magnetic diffusivity drops and the magnetic 
Reynolds number increases above one in a manner very similar to that shown in Figure [S]^ 
(dashed line). The equatorward flows advect field lines from high to low latitudes, so their 
geometry in the vicinity of the base of the convection zone is mostly horizontal. This effect 
is the pre-confinement process discussed in Section [3] 

Interestingly, we find that the magnetic Reynolds number remains relatively high in the 
tachocline, even though the meridional flow velocity drops significantly. Polar field lines 
are bent horizontally, pushed downward and advected toward mid-latitudes (roughly 30°) 



as predicted by GM98 and studied by Wood & Mclntyre (2011). Below the tachopause at 
r = O.5i?0, by contrast, we see that the field is not affected by the flows and is more or less 
equal to the imposed primordial field (see Figure [s]). 

The bottom right panel of Figure |6] shows the angular velocity profile of the reference 
model. The convection zone, as expected, rotates with the imposed differential rotation given 



by Equation (18). Most of the radiation zone (for r < O.Si?©), by contrast, exhibits a uniform 
angular velocity with firz — 0.92r2oq, which extends from the base of the tachocline to the 
core. The fact that fi^z is so close to the observed rotation rate in the solar radiation zone 
is interesting, but could be a coincidence. 

The tachocline lies between these two regions. We note that its rotation rate (in the 
reference model) is not a monotonic function of radius at all latitudes. Indeed, the polar 
region rotates even more slowly than the convection zone above, a feature that is not seen 
in the solar tachocline. This sub-rotating region is caused by the extraction of angular 
momentum by the equatorward meridional flow. Because of our unrealistically high thermal 
diffusivity, the strength of the meridional flow in our model is much larger than expected 
in the Sun, which in turn causes this sub-rotating feature to be unrealistically wide and 
strong. We show in Section [Sjthat extrapolating our results towards more realistic parameters 
leads to weaker meridional flows, and that the amplitude and extent of the sub-rotation 
decreases. In other words, the presence of this feature in our reference simulation should not 
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Figure 7. Comparison of the various terms contributing to the azimuthal vorticity equation in our reference model, in units of 
f2?j. Only the radiation zone is shown. The rotational shear (top left panel) balances the baroclinicity terms (top right panel). 
Both viscous (bottom left) and magnetic (bottom right) vorticity production terms are negligible. 

be viewed as an inherent problem with the GM98 model, but instead, as a natural limitation 
of simulations that have to be run, for computational feasibility, at non-solar parameters. 

The remaining panel in Figure M (top right) shows the temperature perturbation profile 
in the reference model. The latter is negligible in the magnetically dominated part of the 
radiative region, below the tachopause. By contrast, the tachocline exhibits significant radial 
and latitudinal gradients consistent with thermal equilibrium and thermal-wind balance (see 



Section 4.2). 



We now analyse the results of the reference model more quantitatively to determine the 
dominant force balance in each region, and compare our findings with the models of GM98 



and Wood et ah (2011) 



4.2 Thermal wind balance in the radiative interior 

GM98 begin their scaling analysis by assuming that the tachocline and tachopause are in 
thermal-wind balance. We can easily verify this assumption in our model. Thermal-wind 
balance occurs when the system satisfies both hydrostatic and geostrophic equilibrium. Tak- 
ing the curl of the steady-state momentum equation divided by the background density p, 
we obtain: 



2{nQ + n)rsme 
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Figure 8. Comparison of the various terms that contribute to the azimuthal force balance in the radiation zone in our reference 
model, in units of poRqQ.%. In the tachocline and tachopause (for r > O.S-Rq), the Coriolis force (top left) and Lorentz force 
(top right) are in balance, while the non-linear advection term (bottom left) and the viscous torque (bottom right) are negligible. 
Note that the tachocline is mostly force-free while the tachopause stands out as the region where strong Coriolis and Lorentz 
forces balance. 



where the first term measures the rotational shear along the rotation axiqj, and the next 
two terms account for the total baroclinicity. The balance between rotational shear and 
total baroclinicity is the thermal-wind equation. The remaining terms in this equation are 
the magnetic and viscous vorticity-production terms as well as the curl of the gyroscopic 
pumping term. The latter is included for completeness but vanishes in the radiative interior 



(see Section 2.5.1) 



The two top panels in Figure u^ show that the rotational shear (left) and the total baro- 
clinicity term (right) clearly balance each other in the tachocline and tachopause region. The 
bottom panels show the viscous (left) and magnetic (right) terms, and confirm that their 
respective contributions to the azimuthal vorticity equation is negligible. In our solutions, 
thermal- wind balance holds throughout most of the radiative interior, as assumed by GM98. 



In the model of Wood et al. (2011), on the other hand, thermal-wind balance is broken by 



Lorentz forces in the tachopause at the bottom of the tachocline. Our results show that this 
is not necessary in order for a tachopause to form. We note, however, that the strength of 
the primordial magnetic field in our reference model is much weaker than the typical field 



strengths considered by Wood et al. (2011). With a stronger magnetic field, thermal- wind 



balance might well be broken in the tachopause. 



The two terms in the bracket on the left-hand side can be expressed more concisely as e^ • VQ. 
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4.3 The azimuthal force balance in the tachocHne 

We now examine the second ingredient of the GM98 model, namely the angular-momentum 
balance, by inspecting the azimuthal component of the momentum equation. The latter can 
be expanded as: 

— 2pQq {uq cos 6 + Ur sin 6) (Coriolis) 
dUfj, ue dus 1 cos 6 



or r GO r rsmu 



(Inertial) 
+jrBg — jeBr (Lorcutz) 



— [us — Ucz] (Gyroscopic pumping) 

r 

+ (V ■ 11)0 = (Viscous stresses) (31) 

where 11 = pz/[Vu + ( Vu)^ — | ( V ■ u)I] . Because neither gravity nor pressure enters into the 
angular-momentum balance, forces that are insignificant in the meridional directions can be 
of leading order here. Figure [s] shows each of these terms individually (except the pumping 
term, which does not act in the radiation zone). 

First and foremost, note that the viscous stresses are negligible in the tachocline, confirm- 
ing the theoretical expectation that viscosity should play no role in its dynamics [j This is in 
contrast with the results of the direct numerical simulations of Strugarek et al.| (|2011) and 



Rogers (2011), in which a ^ 1 and angular-momentum transport is viscously dominated. 

The top panels of Figure |8] show, from left to right, the azimuthal component of the 
Coriolis and Lorentz forces. Both are in balance in most of the radiation zone (except close 
to the polar axis near the inner boundary), strongest in the tachopause, and significantly 
weaker in the tachocline. This is consistent with the GM98 model, in which the tachocline is 
essentially force-free. Lorentz forces become important only in the tachopause, where they 
provide the angular momentum necessary for the fiows to be defiected equatorward. 

4.4 The thickness of the tachopause and the tachocline 

The force-free nature of the tachocline constrains the fiow of material downwelling from the 



convection zone (see Wood et al. 2011, for a detailed discussion). Indeed, if viscous torques. 



magnetic torques and inertial terms are neglected in Equation ( [3l| ), then 

^ Viscous stresses are significant in the proximity of the polar axis, where they balance the Lorentz forces in a thin Eknian— 
Hartmann-type boundary layer that is produced by the artificial lower boundary at r = rin. Such a boundary layer is not 
present in the Sun. 
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Figure 9. Radial mass flux \pur\ at 80° latitude. A dotted line is used when Ur < 0, and a solid line is used when Ur > 0. 
We observe depth-independent downwelling from the radiative— convective interface down to about r = 0.58i?Q, a region we 
identify as the tachocline. Below O.58-R0, \pur\ decreases to zero then oscillates about it. We identify the region between the 
base of the tachocline and the first zero of \pUr\ as the tachopause (here at r = O.525-R0). 



[Qq X u),^ = 



(32) 



which iniphes that the meridional flow velocity must be parallel to the rotation axis. Com- 
bining this result with mass conservation (in a cylindrical coordinate system) we then have 
d 



dz 



(pu,) = , 



(33) 



where z = rcosO. At high latitudes, Ur — Uz, so the radial mass flux must be roughly 
constant along the polar axis. Within the tachopause, by contrast, the prograde torque from 
the Lorentz force is signiflcant and gyroscopically pumps the fluid equatorward, allowing the 
meridional flow to turn over and return to the convection zone. As a result, pUr is no longer 
constant. 

This qualitative picture is verifled in Figure [ol which shows the mass flux \pUr\ at 80° 
latitude. We see that \pUr\ is constant between the radiative-convective interface and r ^ 
0.58Rq and then decreases rapidly to zero, where the magnetic torque is strongest. 

This observation provides an objective measure of the location of the tachopause, and 
hence of the thickness of the tachocline. In all that follows, we deflne the base of the tachocline 
rx (which is also the top of the tachopause) to be located at the depth where \pUr\ drops 
by 5% from its value at r = Vcz, as measured at 80° latitude. In Figure [9| we flnd that 
tt = 0.575Rq in the reference model. We then deflne the tachopause as the region located 
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Figure 10. The toroidal field in the radiation zone, in units of Bq. Note how the strongest field concentrations are found 
within the tachopause, and within a narrow region near the base of the convection zone around 45° latitude. 

between the base of the tachochne and the depth at which the vertical mass flux first equals 
zero. From Figure |9j we find that this happens at r^ = O.525i?0. The thickness of our 
reference model tachocline is then roughly three times that of the tachopause. 

Although the particular way in which we measure the tachocline and tachopause thick- 
nesses (e.g. the latitude at which the measurement is made, the percentage drop in \pUr\) 
are somewhat arbitrary, the concepts behind the definitions themselves are robust. This 
provides a means to test how the thicknesses of the tachocline and tachopause vary with 
the governing parameters (in particular, E^^ and Ejj) and is used in Section Is] to test the 
predictions of the GM98 model. 



4.5 Toroidal field 

The final aspect of the GM98 theory that we investigate is the distribution of toroidal field in 
the tachocline and in the tachopause. Strong toroidal fields are generated when the poloidal 
component of the magnetic field is twisted by the azimuthal flows. GM98 argue that this 
effect should be most important in the tachopause, which lies at the interface between a 
region of strong angular- velocity shear (the tachocline) and a region dominated by a strong 
poloidal field (the deep interior). They also mention the possibility that toroidal fields could 
be produced in the mid-latitude upwelling region (where the flow drags poloidal field lines 
back into the tachocline), and might thus play a role in the solar dynamo. 



Figure [TO] shows the distribution of toroidal field in the radiative region for our reference 
model. As expected, we find that strong fields are indeed generated in the tachopause. We 
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also find them in a localized region of the tachocline, around 45° latitude, which is somewhat 
closer to the poles than the center of the upwelling region (which lies around 30° latitude). 

5 NUMERICAL EXPERIMENTS VARYING THE PARAMETERS 

Having established that our reference model satisfies the various force balances predicted 
by the GM98 model, we now test the scaling laws implied by these balances. We do so by 
varying the parameters governing the system around the reference model values. 

5.1 The effect of the thermal diffusivity 

We first vary the thermal diffusivity, keeping all other parameters fixed - in particular, 
without changing v or N . We do this by multiplying the reference model function ^(r) (see 



Table 1 and Section 2.3) by a constant factor /«. This implies that K{r) changes in the 
convection zone as well, but it can be shown by inspection of the meridional flow velocity 
profile for r > r^z that this alone has a negligible effect on the solution. Note that varying 
K without varying u or N implies that a changes as well, but in all cases presented we 
made sure that a remains smaller than 1 to guarantee that the effects of viscosity are indeed 
negligible in the tachocline and tachopause. 

We are interested in finding out how varying k in the radiation zone and the tachocline 



affects the dynamics of these regions. The results are shown in Figure 11 The left panel 
shows the effect of reducing k from the reference model by a factor of 2. Because of thermal 
equilibrium, this directly affects the amplitude of the downwelling meridional flows in the 
tachocline and therefore the depth at which these flows confine the field. We see in the 
streamline plot that the tachocline is now thinner than in the reference case (see Figure M 
for comparison). The field lines are also less distorted in the tachocline than in the reference 
model. The right panel shows the effects of an increase in n from the reference model value 
by a factor of 2.25. This time, the meridional flows are stronger, which leads to a thicker 
tachocline, and greater distortion of the magnetic field lines. 

Another consequence of the respectively weaker and stronger meridional flows is a clear 
change in the angular- velocity profile of the polar tachocline. As discussed in Section |4], we 
find that the subrotating region becomes both weaker and smaller as the meridional flow 
velocity in the tachocline decreases. 



We now study these trends more quantitatively. As discussed by Spiegel & Zahn ( 1992 ) 
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Figure 11. Simulations with a thermal diffusivity k that is half that of the reference model (left, /„ = 0.5), and 2.25 times 
that of the reference model (/« = 2.25, right). All other parameter values are otherwise unchanged from those of Table 1. The 
higher diffusivity simulation has a thicker tachocline, and field lines are more strongly distorted by the downwelling flows. 

and GM98, thermal-wind balance and thermal equilibrium imply specific scalings for the 
tachocline. Indeed, thermal-wind balance means that 
dn g df 



2r sin 6^2(7, 



dz rT de 



(34) 



(see Section 4.2) where z = rcos9. If we further assume that the tachocline thickness 
A < i?0 then: 



(35) 



T^ glA ' 

where we have approximated the rotational shear as dCl/dz ^ a^Q/A, where a characterizes 
the level of differential rotation. The latitudinal temperature gradient is approximated as 
dT /do K. IT, where / is a latitudinal wavenumber of order unity. If the tachocline is also in 
thermal equilibrium, then 






g dr 

so that 

^ iV2A2 f 



2 ' 



(36) 



(37) 



30 L. A. Acevedo-Arreguin, P. Garaud & T. S. Wood 

using d/dr ^ /3/A, where /3 is a geometrical constant of order unity. From these two equi- 
hbria, GM98 deduce that 

^^ I N^ re, A3 "" iV2 re. A3 ' ^•'''^ 

where i^ is a geometrical constant. Using a = 0.07, (5 = -k and / = 4.5, they estimate that 
ir~0.15. 



Equation (38 ) relates the amplitude of the downwelling meridional flows in the tachocline 
to its thickness under robust physical assumptions; it is a generic property of many models 
dSpiegel fc Zahnl|1992[ |Gough fc Mclntyre||1998| [Wood et al.||2011| ). We can now check its va- 
lidity using our simulations in conjunction with the measurement technique for the thickness 
of the tachocline and tachopause described earlier. However, the comparison is not entirely 
trivial, since the tachocline and tachopause in our simulations are not necessarily as thin as 
required by the assumptions made by GM98 (recall that our diffusion parameters are still 
significantly larger than those of the Sun). This has two consequences. First, the background 
density and buoyancy frequency vary by a significant amount within the tachocline. This 
implies that although pUr is constant with depth, Uj. is not. One must then choose at which 



radial position Equation (38) should be used to estimate Ur. In what follows, we take for 
simplicity r = 0.7 Rq, which is well within the tachocline for all cases considered. Further- 
more, GM98 assume that the tachopause is much thinner than the tachocline, so they did 
not need to differentiate between A (the tachocline thickness only) and A + 6 (the sum of 



the tachocline and tachopause thicknesses). Their original derivation leads to Equation (38). 
Here, 6 can be of the same order of magnitude as A, so we must decide whether Equation 



(138]) or 

02 3 

Ur ^ Kut^eor = K^ ^ ^ whcrc D = A + S (39) 

is more appropriate. We actually believe that the latter is, since the tachopause is also 
in thermal equilibrium and in thermal-wind balance, and since both the shear and the 
temperature perturbations should vanish at the base of the tachopause rather than the base 
of the tachocline (so dVt/dz ~ uVLq/D and dT/dr ^ (3T / D instead). 

The relationship between Z^, the amplitude of the mass flux in the tachocline \pUj.\, and 



the thickness of the tachocline-tachopause region D = 5 + A is shown in Figure 12 As 
expected, increasing /^ effectively weakens the stable thermal stratification, allowing for a 
larger meridional mass flux in the tachocline. We find that pUr changes by a factor of about 
two when /^ changes by a factor of five. The stronger flows push the field down further so D 
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Figure 12. Variation with respect to /« of the mass flux \pur\ in the tachocHne (dotted line), of the combined thicknesses of 
the tachochne and tachopause D (dashed hne), and of ii" = |«r|/«thcor (solid line), all evaluated at r = 0.7Rq. Our results 
show that Equation ||39|l provides a good estimate for |ur|, and that K ~ 0.5. 



increases too, although by a smaller amount (about 25% only). Using this information, we 



can check the estimate for Ur given in Equation (39). This is also shown in Figure 12 We see 



that the ratio K = |ur|/Mtheor varies very little with f^, taking a value of about 0.5. This is a 
few times larger than the estimate provided by GM98, but is overall remarkably consistent 



with their predictions. Finally, note that if we use A instead of D, that is. Equation (38) 



instead of (39), the ratio |Mr|/wtheor varies more with /^j. This suggests that Equation (39) 
is the better estimate for Ur- 



In summary, we find that our simulations confirm the estimates of GM98 and Wood et al 



(2011) concerning the relationship between the downwelling flow velocities and the thickness 



of the tachocline, albeit with a minor modification. However, we note that the range of /„ 
tested is very limited, owing to computational limitations. With the other parameter values 
and resolution fixed, we were unable to find any solution for larger /^ than 2.25, and for 
smaller /^ than 0.5. For larger f^, the tachocline is increasingly deep and the meridional 
flows interact more strongly with the magnetic field. Our steady state solver has difficul- 
ties following the solutions in that limit. For smaller /k, a increases above unity and the 
tachopause begins to overlap with the convection zone, both of which lead to a qualitative 
change in the nature of the solutions. 
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5.2 The effect of the magnetic diffusivity 

We now study the effect of varying the magnetic diffusivity of the radiation zone on the 
dynamics of the system. Starting from the reference model described in Section [4], we vary 
r/rz (or equivalently, E^), by multiplying the reference value (see Table 1) by a factor /^. 
All other parameters remain the same as in the reference model. This time, the magnetic 
diffusivity profile within the convection zone does not change when varying /^, which implies 
that the magnetic Reynolds number within the convection zone remains approximately the 
same across all simulations. This ensures that the pre-confinement process is the same for 
all cases. 



Figures 13 1 and 13 d show our results for /^ = 0.5 and /^ = 3.33 respectively. The 
low- diffusivity case clearly exhibits a much stronger degree of field confinement within the 
radiation zone than the high-diffusivity case (see lower left panels). From the streamline 
plots (top left panels), we also see that the tachocline is slightly thicker when /^ is lower. 
Figure [14^ shows that the geometry of the tachocline and the tachopause, estimated using 



the method described in Section 4.4, changes significantly with /^. As /^ decreases, we 
find that the base of the tachocline moves deeper into the radiation zone, but the base 
of the tachopause remains roughly at the same place. This implies that the tachopause 
becomes thinner whereas the tachocline thickens, and the ratio 5/A decreases. At low enough 
magnetic diffusivity (e.g. for /,, = 1 and lower), the numerical solution satisfies the condition 
5 < A assumed by GM98. 

A defining property of the tachopause is balance between advection and diffusion of the 
magnetic field. GM98 argue that, in a steady state, 

\Ur\ ^ ^ , (40) 

which relates the downwelling meridional flow velocity to the thickness of the tachopause 
and the magnetic diffusivity. In other words, the magnetic Reynolds number within the 
tachopause, Rm^ = 5\ur{rT)\/rj-[z, where UrivT) is the radial velocity at the base of the 
tachocline, should be of order unity. To study whether this scaling holds here, we compute 



Rjqt as a function of /^. The results are plotted in Figure 14 d. We see that, although Ur 
itself changes by a factor of about three, Rm^^ is essentially constant and close to unity. This, 
once again, is in good agreement with the GM98 model. Finally, we see that, as a natural 
consequence of the decrease in Ur when /^ decreases, the strength and extent of the polar 
sub-rotating region also decreases, as discussed in Section |4j 
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Figure 13. Simulations with a magnetic diflfusivity rjrz that is half that of the reference model (left, /,, = 0.5), and 3.33 times 
that of the reference model (fn = 3.33, right). All other parameter values are otherwise unchanged from those of Table 1. The 
lower-diffusivity simulation has a thicker tachocline, and the magnetic field lines are more distorted. 
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Figure 14. Left; Variation of the position of the base of the tachocline (dotted line with symbols) and tachopause (dashed line) 
with magnetic diffusivity, as measured by /^. Right: While the downwelling flow velocity steadily increases with /^ (dashed 
line), the magnetic Reynolds number of the tachopause, defined as Rm-r = |«r(»"T)|<5/^rz, remains close to unity (solid line 
with symbols). 
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5.3 The influence of the Elsasser number 

The results discussed so far prove that the assumptions exphcitly made by GM98 are correct, 
and strongly suggests that their model is a good representation of a laminar solar tachocline. 
The scalings we have been able to test appear to hold, at least within factors of order unity. 
Within the scope of that model, GM98 deduce that the thickness of the tachocline should 
vary with the primordial magnetic field strength as: 

where the constant of proportionality between the two sides of this expression also depends 
on a, /3 and /. 

Unfortunately, we have not been able to test this last relationship directly. The weak 
dependence on the internal field strength implies that one would need to vary Bq (or equiv- 
alently. A) by many orders of magnitude to detect any significant change in A. Here, with 
all other parameters fixed at their reference model value, we have only been able to increase 
Bq by about one order of magnitude, and, within the range of Bq tested, have not detected 
any significant variation in A. Any attempt to increase Bq further results in the failure of 
our relaxation algorithm to find steady-state solutions. We believe the source of the problem 
is similar to the one encountered by GG08; as Bq increases, the system probably becomes 
unstable to various MHD instabilities associated with unrealistically large field strengths 
near the artificial inner boundary of the domain, and a steady-state solution can no longer 
be found. Numerically speaking, the matrix inversion problem involved in the steady-state 
search becomes increasingly ill-conditioned, and the solver no longer converges. By selecting 
other values of the diffusivities, it is possible to find solutions for larger field amplitudes. 
However, in no cases have we been able to span a range for Bq large enough to test the 



validity of Equation (41), for the same reasons as the ones described above. Nevertheless, 
since Equation (41) is possibly the least robust result of the GM98 model anyway (see Sec- 
tion [6] for details), we do not view our failure to validate it directly in our simulations as 
particularly troublesome. 

Beyond interfering with the relaxation algorithm in our simulations, one should cer- 
tainly question whether hydrodynamic or magnetohydrodynamic instabilities could actually 



prevent the GM98 model from applying to the Sun altogether. Brun & Zahn (2006) inves- 
tigated this issue, with their 3D time-dependent model of the radiation zone. They found 
two emerging modes of instability. MHD instabilities with high azimuthal wavenumber m 
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seem to develop deep in the radiation zone, and similar features are reported by Strugarek 



et al. (2011). As discussed by both studies, these are in fact expected since the simulations 



are initiated with a field in a purely poloidal configuration, which is known to be unstable 



to such modes (Wright 1973 Markey & Tayler 1973). Brun & Zahn (2006) also report that 



instabilities with a characteristic m = 1 azimuthal wavenumber (Tayler 1973) are found in 
regions with large toroidal fields, in particular close to the pole. 



Interestingly, however, both Brun & Zahn (2006) and Strugarek et al. (2011) find that the 



deeply-seated, high-m instabilities do not disrupt significantly the structure of the internal 
poloidal field. The growing instabilities generate a non-axisymmetric toroidal field whose 
addition rapidly stabilizes the system ( Braithwaite fc Spruit|2004 ), and there is no indication 
that other MHD instabilities play any further significant role in the dynamics. From this we 
conclude that even though they could be the reason behind the failure of our steady-state 
solver to find solutions for larger field strengths and lower diffusivities, MHD instabilities do 
not necessarily imply a failure of the GM98 model itself. The potential effect of turbulence 
and other instabilities on the GM98 model is discussed in more detail in Section |6l 



5.4 The effect of the viscosity and the failure of previous models. 

As described in Section [T| attempts at finding numerical representations of the GM98 model 
using large-scale 2D or 3D time-dependent simulations have been made by |Brun &: Zahn 



(2006), Strugarek et al. (2011) and Rogers (2011). In all three cases, the attempts failed to 



exhibit any evidence for confinement of the magnetic field below the base of the convection 
zone. 



Strugarek et al. (2011) concluded that the 2D laminar nature of the GM98 model is 



too simplistic, and cannot capture the dynamics of the solar tachocline adequately. 

We believe this conclusion is premature. While the GM98 model neglects any effects 
of turbulence and instabilities (see Section |6] for detail), its 2D laminar nature is not the 
reason for the failure of previous models. A strong clue to this effect comes from the analysis 



performed by Strugarek et al. (2011) and Rogers (2011) of their own simulations, which 



demonstrates that momentum transport therein is essentially viscously dominated, whereas 
viscosity is assumed to be negligible in the GM98 model. These numerical simulations clearly 
operate in a very different physical regime from that envisaged by GM98. 

We argued earlier that the source of the problem lies in the parameter a, which is ^ 1 



in the work of Strugarek et al. (2011) and Rogers (2011) but < 1 in the solar tachocline. To 
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Figure 15. Left: DifEusivity profiles actually used by |Strugarek et al.| | |20rT| |, as well as one where r;(r) is shifted upward to 
allow pre-confinement (see main text for detail). Right: The a profile used by |Strugarek et al.| l [2011[ |. Note how it is 3> 1 
everywhere in the radiation zone. 



demonstrate this effect more directly, we now present an axisymmetric steady-state solution 



using very similar parameters to those of Strugarek et al. (2011), and compare our results 
with theirs. 



Following Strugarek et al. (2011), we use solar profiles for all background thermody- 



namical quantities, including the buoyancy frequency N{r). We also use their diffusivity 
profiles: 

V-O.6753i?0~ " 



z/(r) = 8.0 X 10^ + 8.0 x 10^^ 
r^{r) = 8.0 X 10^° + 1.6 X 10^^ 
K(r) = 8.0 X 10^2 + 3.2 x 10^^ 



1.0 + tanh 
1.0 + tanh 
1.0 + tanh 



O.Oli?0 
r - O.6753i?0 

OMRq 
r - O.6753i?0 

O.OlRr. 



(42) 



which are given in cgs units here, and shown in non-dimensional form in Figure 15 a,. Note 



that as a result of this choice, the corresponding profile for crsBziil^); shown in Figure 15 3, 
is quite different from o"0(r). In particular, o"sbzii ^ 1 in the entire radiative region, while 
(70 should be smaller than one in the tachocline. 



To emulate the 3D, time-dependent simulations of Strugarek et al. (2011) with our 



steady-state solver, we use a similar method to the one presented in Section |2j We impose the 
differential rotation in the convection zone using the same forcing term as in our reference 



model (see Section 2.5.1). We must also assume the presence of an azimuthal current to 



maintain the primordial magnetic field, as discussed in Section 2.5.2 We take its functional 
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form to be the same as in our reference model, and choosq^ A = 6 x 10^. Note that since 
the magnetic field continuously decays by Ohmic diffusion in a time-dependent calculation, 
it is difficult to know exactly what value of Bq most appropriately describes the simulation 



of Strugarek et al. (2011) at any point in time. As such, the following comparison between 



their results and ours necessarily remains mostly qualitative. 

The resulting steady-state solution for the system described above is shown in Figure 
16 a.. It is both strikingly different from the solutions presented in Section |4| and quite similar 



to the final stages of evolution of the simulation presented by Strugarek et al. (2011). The 



downwelling flows entering the radiation zone from above decay rapidly within the interior 
and are clearly unable to confine the internal magnetic field. The latter diffuses into the 



convection zone, and, by Ferarro's isorotation law (Ferraro 1937), causes the radiation zone 



to rotate differentially. We also see that the meridional flow pattern just below the base of 
the convection zone is very different from that of our reference model (see Figure p«, but 



rather similar to that found by Strugarek et al. (2011) and by Rogers (2011), with shallow 



alternating layers of opposite circulation. Hence, despite being qualitative, the comparison 



between our 2D steady-state model and the 3D time-dependent model of Strugarek et al. 



(2011) appears to be meaningful. 



Of course, one may correctly argue that our reference model facilitates the field confine- 
ment process through jore-confinement (see Section [3]). To test whether pre-confinement can 
produce a confined magnetic field even for a ^ 1, we present a second simulation, set up 
exactly as above, with the exception of the magnetic diffusivity profile. The latter is now 



shifted upward, as shown in Figure 15 1, thus allowing for pre-confinement to take place. The 



resulting 2D steady-state numerical solution is shown in 16 d. Even though the magnetic field 
is now indeed somewhat pre-confined in the convection zone, the field lines in the radiation 
zone remain mostly unaffected by the downwelling fiows, and significant differential rotation 
persists. In other words, pre-confinement is not sufficient, on its own, to produce a confined 
magnetic field at these parameter values. 

We now analyze the results of this second simulation, which includes pre-confinement, 



more quantitatively. Figure 17 shows the relative amplitude of the various forces operating 
in the azimuthal component of the momentum equation. It is clear that, by contrast with 



^ This large value is needed in order to have an Elsasser number that is comparable to that of [Strugarek et al.| | |201l| in the 
region of the tachocline. Indeed, A is defined according to the magnetic field strength near the core, while the amplitude of the 
primordial field rapidly drops with radius so that its value in the tachocline region is much smaller than Bq (see Appendix A). 
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Figure 16. Steady-state solutions obtained in a parameter regime close to the one used by |Strugarek et aI.| ( |20lH , without 
(left) and with (right) pre-confinement (see main text for detail). 



our reference model, the viscous stresses are dominant throughout the radiation zone. They 
balance the Coriolis force at high latitudes, and the Lorentz force at mid- to low latitudes. 



This is consistent with the results of Strugarek et al. (2011), who observe viscous stresses to 



be significant in their simulations, even though the Ekman number is <^ 1. 

In summary, our axisymmetric laminar and steady-state models of the solar interior, 
yield solutions that are comparable to the long-term evolution of more realistic 3D time- 
dependent simulations, when run at the same parameters. Our results demonstrate that 
the failure of prior numerical attempts to model magnetic confinement and a self-consistent 
tachocline can be attributed to the fact that cr ^ 1 in all those simulations. In this "high 
a''^ regime, viscosity plays a significant role in the transport of angular momentum, the 
meridional flow pattern consists of alternating shallow cells with a flow velocity that decays 
exponentially with depth, and the magnetic field is unconfined. 



Dynamics of the solar tachodine 39 




5.10"'' 



5.10" 



0.2 0.4 0.6 0.6 0.2 0.4 0.6 0.8 
r/Ro --/Ro 



Figure 17. Comparison of t he v arious terms that contribute to the azimuthal force balance in the radiation zone for the 
simulation presented in Figure Puip, in units of pq-Rq^!?,. Shown are the Coriolis force (top left), Lorentz force (top right), the 
inertial term (bottom left) and the viscous torque (bottom right). By contrast with our reference model, the dominant balance 
below the base of the convection zone is between the Lorentz force and the viscous torque. 



6 DISCUSSION AND PROSPECTS 

6.1 Summary 

The GM98 model provides a plausible explanation for the uniform rotation of the bulk of 
the solar radiation zone, and various other observations related to the tachocline, including 
its angular velocity profile, sound speed anomaly, and the rate of lithium depletion. It also 
provides scaling laws that relate the thickness of the tachocline and its mixing timescale to 
the strength of the assumed internal primordial field. However, the GM98 model makes a 
number of simplifying assumptions, and describes only the laminar, axisymmetric, and time- 
independent dynamics of the tachocline. As such, its predictions can only be fully verified by 
self-consistent, nonlinear numerical modeling. Unfortunately, recent attempts at simulating 
the solar interior have failed to reproduce the dynamical regime predicted by the GM98 



model, which have led some to conclude that it is unworkable (Strugarek et al. 2011) . 



In this work, we also study the GM98 model numerically, still assuming an axisymmetric 
steady state, but solving the full set of governing nonlinear MHD equations in global spherical 
geometry. Although our model cannot self-consistently describe all the effects of instabilities 
and turbulence in the tachocline, it can certainly be used to validate (or invalidate) the GM98 
model under these more general conditions. The limitations of our approach are discussed 
in depth below. Within this framework, however, we not only present the first numerical 
simulations of the solar interior to exhibit a tachocline and a tachopause with properties 
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that are qualitatively and quantitatively consistent with most aspects of the GM98 model, 
but also provide a simple explanation for the failure of previous models to do so. 



For low-enough magnetic diffusivity (see Figures ^ and 13 1 for instance), we obtain nu- 
merical solutions of the governing equations that have a radiation zone held close to uniform 
rotation by a confined primordial magnetic field. As in the GM98 model, we find that the 
field is confined by meridional flows downwelling from the convection zone, within a thin 



tachopause in magnetic advection-diffusion balance (i.e. with Rm ~ 1, see Figure 14). The 
tachopause is also in magnetostrophic balance (see Figure Isl), and the Lorentz forces gener- 
ated as the field is wound up by the rotational shear provide just enough angular momentum 
to deflect to the downwelling flows equatorward on their path back to the convection zone. 
As assumed by GM98, the bulk of the tachocline is "magnetic free" (in the sense that the 
magnetic forces are negligible) and in thermal- wind balance (see Figure ItI). The strength of 
the downwelling flow is determined by this balance together with thermal equilibrium (see 



Figure 12). 



By contrast, when we use the same governing parameters as Strugarek et al. (2011 ), we do 



not obtain a confined magnetic field, or a tachocline. As discussed in Section [T| the problem 



is linked to the parameter a = vPrA^/fi©, which is ^ 1 in the models of Strugarek et al 



(2011 ) (as well as Rogers (2011 ) and Brun & Zahn ( 2006[ )), but is < 1 in the solar tachocline. 



When 0" ^ 1, meridional flows are strongly suppressed beneath the convection zone, and 
therefore unable to confine the field, and the angular-momentum balance is dominated by 
viscosity. When a < 1, however, meridional flows are able to extend (or "burrow") much 
more deeply into the radiation zone, and with an amplitude sufficient to confine the magnetic 
field. 

Within the scope of our axisymmetric, steady-state model, our results confirm the pre- 
dictions of GM98, and guide parameter selection for future 3D numerical simulations. We 
now discuss the caveats and implications of this model and of our findings in more detail, 
and lay out future theoretical and observational prospects for solar and stellar astrophysics. 



6.2 Model caveats and their implications 

The axisymmetric, steady-state approach to studying the solar tachocline proposed by GM98 
and adopted throughout this paper has clear numerical advantages, but neglects by construc- 
tion any time- dependent, non- axisymmetric dynamics. In what follows, we loosely refer to 
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these neglected dynamics as "turbulence", although they may also include any effects of 
Alfven waves or internal waves. We now consider whether turbulence in and around the 
tachocline can qualitatively affect the scaling laws derived by GM98. 

6.2.1 The validity of thermal-wind balance and thermal equilibrium 



The main tachocline scaling law given in Equation (39) results from the assumptions of 
thermal-wind balance and thermal equilibrium. As we now show, this scaling probably con- 
tinues to hold regardless of the presence or absence of turbulence, given reasonable upper 
bounds for the turbulent intensity. 

Indeed, one can estimate under which conditions a turbulent flow with energy-bearing 
eddies of vertical scale Ir and horizontal scale Ih can be present in the tachocline without 
upsetting either balance. First note that in strongly stratified flows Ir is typically much 
smaller than //j. If Wj-ms is the mean eddy velocity, then their typical horizontal velocity will 
be of the order of Urms whereas the vertical velocity Wrms is of the order of Wrms ~ {lr/lh)urms- 
Second, note that the momentum (or vorticity) equation involves the spatial derivative of the 
Reynolds stress tensor. This derivative involves vertical scales of the order of the thickness 
of the tachocline A, and horizontal scales of the order of the solar radius Rq. Using these 
estimates, it can then be shown that for thermal-wind balance to hold, one must have 



IV X (V ■ (uu))|^ ~ ^ max [J-, ^j « 2Q^R^^^ ~ IQ-'h-' . (43) 

Finally, note that the development of any instability in the strongly stratified tachocline 
is conditional on the fact that their vertical scale be small enough not to be constrained by 
the buoyancy restoring force. In other words, vertical turbulent fluid motions are by nature 
thermally diffusive, and do not contribute much to the vertical heat fluxj^ The thermal 
energy Equation ^ is thus expected to hold in the presence of self-consistently generated 



turbulence. Using A ~ 0.02/?© in (43), we then find that Equation (39) holds as long as 
«rms < /, . X cm^s^ . (44) 

Meanwhile, a strict upper limit on the turbulence intensity Uj-^s comes from assuming 
that it cannot be larger than a fraction of the convective velocities in the overshoot layer, 
so that Mrms ^ Ncz{0.1Hp) ~ 7 X lO^cm/s, where Ncz ~ lO^^s^"*^ is the imaginary part of the 
buoyancy frequency just above the base of the convection zone, and where we have assumed 

® Helioseismic observations also confirm this statement. 
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that the overshoot layer height is no larger than Q.lHp ~ O.OI-R0, where Hp is the pressure 
scaleheight at the same location. In other words, we expect Wrms to be at most a few meters 
per second, and to decrease rapidly with depth below the base of the overshoot layer. This 



implies that (44) is readily satisfied, and that Equation (39) can be used with confidence to 



model the dynamics of the tachocline. 

6.2.2 The effect of turbulence on the angular momentum balance in the tachocline and 
tachopause. 

By contrast with thermal-wind balance, angular-momentum balance is much more easily 
disrupted by turbulent motions since the only source of angular-momentum transport con- 
sidered by GM98 and in this paper is advection by slow, large-scale meridional flows. In 



many tachocline models, such as the one proposed by Spiegel & Zahn (1992), turbulence 
is assumed to dominate the transport of angular momentum throughout the tachocline. If 
the turbulent transport is very efficient, and also frictional (i.e. down-gradient in angular 
velocity), then it could by itself maintain a thin tachocline even in the absence of a pri- 
mordial magnetic fielcr] In that case the scaling law (39) would still hold, but any scalings 



constructed from angular-momentum balance, in particular Equation (41), would be invalid 



However, for reasons discussed by GM98 and Mclntyre (2007[) (see also Tobias et al. 



2007), turbulence in the tachocline is unlikely to be frictional. As a result, even if it con- 
tributes significantly to the transport of angular momentum in the tachocline and in the 
tachopause, the thickness of the tachocline must ultimately still be constrained by the 



strength of the magnetic field. In that case. Equation (41) should be replaced by a law 
that incorporates the contribution from turbulence to the angular-momentum balance, but 
the qualitative tachocline picture would be unchanged. 

6.3 Theoretical and observational prospects 

6. 3. 1 Computational prospects 

As discussed above, a complete validation of the GM98 model can only come from 3D 
simulations. Clearly, 3D numerical models are still far from being able to achieve solar 

^ We note, however, that turbulent transport in the tachocUne cannot explain the uniform rotation of the rest of the radiation 
zone while the Sun is undergoing spin-down. This by itself provides strong evidence for the existence of a primordial magnetic 
field (iMestel & Weiss|1987| GM98). 
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Figure 18. A 2D slice of parameter space showing the variation of the Prandtl number and N"^ /Q,'^ in the region of the 
tachocline (i.e. for r £ [0.6, 0.7]-Rq), in our reference model simulation, in the model of |Strugarek et al.| | |201l| l and in the Sun. 
The oblique line represents the limit VrN^ /Q? = 1. Lines of constant a are parallel to this line. Simulations above this line are 
viscously dominated, while simulations below this line are not, and thus reside in the same region of parameter space as the 
Sun. 

parameter values. The question thus arises of whether it is possible to model solar-like 
dynamics with non-solar parameters. In this work, we have shown that it is possible, as long 
as model parameters are chosen such that a = VFiN /Vl < 1 everywhere in the tachocline. 



Figure [18] illustrates our proposed route to ensure that this condition is always satisfied. 
While our reference model is not yet at solar parameters, it is at least in the correct "solar" 
region of parameter space. As computational power continues to increase, it will become 
possible to decrease the Prandtl number further. If one gradually increases N back towards 
its original solar value at the same time, thus ensuring that a remains close to its solar value, 
then every simulation along the path described by the green arrow remains in the correct 
region of parameter space, and numerically-motivated asymptotic scalings (to confirm or 



replace Equation (41 ) for instance) could in principle be derived. We shall present the results 



of 3D simulations of the full Sun, following the path described above, in the future. 



6.3.2 Observational prospects 

Meanwhile, the partial validation of the GM98 model presented here opens interesting 
prospects for applications of this theory to other solar-type stars. In particular, it can be 
used to deduce a number of conditions under which a tachocline could be expected to exist, 
and scaling laws for mixing of chemical species and angular-momentum transport within. 
First and foremost, the star must host a primordial magnetic field of sufficient ampli- 
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tude to impose uniform rotation. For this to be the case, the field must be strong enough 

everywhere in the radiation zone to satisfy A(r) > 1, that is, 

B'^{r,e)>A-K7]{r)p{r)Vt^ (45) 

where VL^, is the star's mean rotation rate, and ri{r) and p{r) are its magnetic diffusivity 
and density profiles. The magnetic field plays no significant role in the long-term dynamics 
of stars that do not satisfy this constraint, as found by Rogers (2011) (see Section [I]). 



Whether (45) holds or not at any point throughout stellar evolution depends primarily on 
two competing processes. On the one hand, magnetic braking spins the star down, and 
weakens the constraint over time. On the other hand, the internal field also decays with 
time through Ohmic dissipation. 

While the rotation rate of a star can be measured, its internal field strength cannot, and 
one must rely on rough arguments to estimate its amplitude. Since the Ohmic dissipation 



timescale for a confined dipolar field is of the ordei^° of tohm = ^Iz/^'^^i where f/ is a mean 
value of the magnetic diffusivity across the radiation zone, and ^ some geometrical factor, a 
necessary condition for the field to be non-negligible is tohm > t* where t^ is the age of the 
star. 

Next, one must study whether the condition a < 1 holds or not. If it does, then the field 
is likely to be "deeply" confined by the large-scale meridional fiows, as discussed by GM98. 
If it does not, the field may still be confined, but this time only by turbulent magnetic 
pumping. The implications for the depth of the mixed layer and its ventilation timescale in 
both cases are quite different. If a > 1, the mixed layer is likely as deep as the convective 
overshoot layer - but no deeper - and its mixing timescale is very short. This implies that 
angular-momentum transport between the convection zone and the radiative interior is very 
efficient (with implications for the stellar spin-down problem), and that chemical species are 
efficiently mixed down to the base of the overshoot layer, but probably no deeper. 

If 0" < 1 on the other hand, the star is likely to have a tachocline that has a similar 
structure to the one described by GM98 and in this paper, albeit probably with a differ- 
ent tachopause (see the previous Section). Since the (unknown) details of the tachopause 
structure relate the thickness of the tachocline to the (unknown) internal field strength, it 

^^ It can be shown that for a dipolar field confined in a sphere of radius rcz which has a uniform magnetic diffusivity rj, then 
iohm = ^cz/?^*? where ^ ~ 4.49 is the first zero of the J3 /2 Bessel function. For variable diffusivity, a good approximation for 
f; is f; = (/J^"^ r]^i'^{T)dr^ (Gough, personal communication, see also 



Garaudll999 1 
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is difficult to place any constraints purely from theory on the plausible depth of a stellar 
tachocline. However, if the latter can be measured via asteroseismology, then the velocity 



of the large-scale meridional flows mixing it can be estimated using Equation (39). This 
estimate should then be used, together with the measured mixed-layer depth, to infer the 
expected surface abundance of various elements. Comparing these predictions with obser- 
vations can then serve as an independent validation (or invalidation) of the GM98 model. 
Finally, the associated estimate for the angular-momentum transport timescale across the 
tachocline can be used to study stellar spin-down (see Oglethorpe & Garaud, in prep.), 
and statistical comparison with observations of the angular-velocity distributions of various 
stellar clusters can serve as yet another independent test of the theory, as in the work of 



MacGregor & Brenner (1991) and Denissenkov (2010) for instance. 



To conclude, our results lay the foundations for future quantitative comparisons between 
models and observations, of the differential rotation profile of course, but also of the light- 
element abundances and of the mult i- dimensional sound-speed structure of the solar interior. 
It will also provide better dynamical constraints, and perhaps new paradigms, for solar 
dynamo theory. Finally, it lays a clear path forward in the study of the internal dynamics 
of other solar-type stars, with natural implications for asteroseismic data from CoRoT and 
KEPLER. 
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APPENDIX A: DERIVATION OF THE BACKGROUND MAGNETIC 
FIELD PROPERTIES 

The expression for the primordial magnetic field Bq generated by the imposed azimuthal 



current jo is most easily derived by introducing the magnetic potential Ao{r, 6) (see Garaud 



k Guervilly 2009) defined as 
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which can be solved to yield 

if ^ r ^ Tq , 

[js^^ - ^{ra + ny + \ranr + f + ^] sin'^ e if r„ ^ r ^ n , (A5) 
sin^ 9 ii rt, ^ r . 

Note that some of the integration constants were set to zero to guarantee that the field in 
the core and at infinity be regular. This then implies that Bq = BqBz for r < r^, where Bq 
is related to Jq (see below). 

The constants Jq, ci, C2, and C3 are obtained by requiring continuity of Aq and dAo/dr 
at r = Tq and r = r^, which results in: 
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Note that this last expression can also be derived directly and more easily through the use 
of Biot-Savart's law. 



We can then derive the components of the background magnetic field Bq using (Al ). We 
find, for instance, that its radial component is 

Bq cos 6 if ^ r ^ r^ , 



Bor{r,i 



'^-W'ils'''' -U''- + ''by + l''-''br+f + ^] cos 9 iira^r<:n, (AlO) 
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if rj, ^ r . 
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Note that the solution guarantees that the amphtude of the magnetic field at the center of 
the sphere to be finite, an advantage from a computational perspective over the expression 
for a point-dipole used by GG08. 
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